Ground state approximation for strongly interacting systems in arbitrary dimension 
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We introduce a variational method for the approximation of ground states of strongly interacting 
spin systems in arbitrary geometries and spatial dimensions. The approach is based on weighted 
graph states and superpositions thereof. These states allow for the efficient computation of all 
local observables (e. g. energy) and include states with diverging correlation length and unbounded 
multi-particle entanglement. As a demonstration we apply our approach to the Ising model on ID, 
2D and 3D square-lattices. We also present generalizations to higher spins and continuous-variable 
systems, which allows for the investigation of lattice field theories. 
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Strongly correlated quantum systems are of central in- 
terest in several areas of physics. Exotic materials such 
as high- Tc superconductors and quantum magnets ex- 
hibit their remarkable properties due to strong quantum 
correlations, and experimental breakthroughs with e.g. 
atomic gases in optical lattices provide a perfect play- 
ground for probing strongly correlated quantum systems. 
The main obstacle in understanding the behavior of those 
quantum systems is the difficulty in simulating the ef- 
fective Hamiltonians that describe their properties. In 
most cases, the strong correlations in the exponentially 
large Hilbert space render an exact solution infeasible, 
and attacking the problem by numerical means requires 
sophisticated techniques such as quantum Monte Carlo 
(QMC) methods or the density matrix renormalization 
group (DMRG) approach 0,0]. 

QMC methods suffer from the sign problem which 
makes them inappropriate for the description of fermionic 
and frustrated quantum systems. DMRG is a variational 
approach that provides approximations to ground states, 
thermal states and dynamics of many-body systems. Re- 
cent insight from entanglement theory have lead to an 
improved understanding of both the success and the lim- 
itations of this approach. Indeed, the accuracy of the 
method is closely linked to the amount of entanglement in 
the approximated states |^,^. Matrix product states [3, 
which provide the structure underlying DMRG, are es- 
sentially one-dimensional and the entanglement entropy 
of these states is limited by the dimension D of the ma- 
trices, which in turn is directly linked to the computa- 
tional cost 0, . Hence a successful treatment of sys- 
tems with bounded entanglement, e.g. one-dimensional, 
non-critical spin systems with short range interactions, 
is possible, while the method is inefficient for systems 
with an unbounded amount of entanglement, e.g. critical 
systems and systems in two or more dimensions. Promis- 
ing generalizations that can deal with higher dimensional 
systems have been reported recently However, the 



computational effort and complexity increases with the 
dimension of the system. In addition, the amount of 
block-wise entanglement of the states used in Ref. 0| 
still scales proportional at most to the surface of a block 
of spins, whereas in general a scaling in proportion to 
the volume of the block is possible. Such a scaling can 
in fact be observed for disordered systems j3j or systems 
with long-range interactions 9]. 

Here we introduce a new variational method using 
states with intrinsic long-range entanglement and no bias 
towards a geometry to overcome these limitations. We 
first illustrate our methods for spin-1/2 systems, and 
then generalize them to arbitrary spins and infinite di- 
mensional systems such as harmonic oscillators. In finite 
dimensions, the method is based on a certain class of 
multiparticle-entangled spin states, weighted graph states 
(WGS) and superpositions thereof. WGS are a 0{N^) 
parameter family of A^-spin states with the following 
properties: (i) they form an (overcomplete) basis, i.e. 
any state can be represented as a superpositions of WGS; 
(ii) one can efficiently calculate expectation values of any 
localized observable, including energy, for any WGS; (iii) 
they correspond to weighted graphs which are indepen- 
dent of the geometry and hence adaptable to arbitrary 
geometries and spatial dimensions; (iv) the amount of en- 
tanglement contained in WGS may be arbitrarily high, 
in the sense that the entanglement between any block of 
iV^ particles and the remaining system may be 0{Na) 
and the correlation length may diverge. 

Note that (iii) and (iv) are key properties in which this 
approach differs from DMRG and its generalizations and 
which suggest a potential for enhanced performance at 
least in certain situations, while (ii) is necessary to effi- 
ciently perform variations over this family. In the follow- 
ing we will outline how we use superpositions of a small 
number of WGS as variational ansatz states to find ap- 
proximations to ground states of strongly interacting spin 
systems in arbitrary spatial dimension. 
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Properties of WGS. WGS are defined as states of N 
spin-1/2 (or qubits), that result from applying phase 
gates Uabifab) — diag(l, 1, 1, e^*'^"'') onto each pair of 
qubits a,b £ {l,2,...,iV} of a tensor product of a^- 
eigenstates |+) = (|0) + |l))/\/2, followed by a single- 
qubit filtering operation Da = diag(l, e'*"), da G C and a 
general unitary operation Ua 

N N 

\^r,d,u)(xl[UaDa n Uati^ab)\+)^^ ■ (1) 
a=l b=a+l 

The phases ipab can be associated with a weighted graph 
with a real symmetric adjacency matrix Tab — 'fab- 
For convenience, we define a deformation vector d = 
(di, d2, ■ ■ ■ , d]\r) and U = (^^C^q. The deformations make 
WGS as used in this letter slightly more general than the 
WGS used in Refs. where da — 0. One can conve- 

niently rewrite |^'r,d,c/) as 

s 

where the sum runs over all computational basis 
states, which are labelled with the binary vector s = 
(si, S2, . . . , Sat)"'". Our class of variation states comprises 
superpositions of WGS of the form 

m 

l*> «E"»l*r,dW,y), (3) 

i=l 

i. e. the superposed states differ only in their deformation 
vector d^*' , while the adjacency matrix T and the unitary 
U are fixed. Such a state is specified by N{N — l)/2 -|- 
37V -I- 2(iV + l)m = 0{N'^) real parameters. 

We now proceed to verify the properties set out in the 
introduction. For property (i), observe that for any fixed 
r and U, all possible combinations of Da G {ai°'\l^'^^} 
lead to an orthonormal basis (note that com- 
mute with Uab)- Hence any state |^') can be written in 
the form Eq. for sufficiently large m < 2^, which 
shows the exhaustiveness of the description. 

The relevance of employing deformations lies in the ob- 
servation that only of the form of Eq. (jSJ permit the 
efficient evaluation of the expectation values of localized 
observables A, i.e. satisfy property (ii). For simplicity 
we restrict our attention to observables of the form 

A = Y,Aab + Y.^-' (4) 

a<b a 

where Aab has support on the two spins a, b. The 
method can be easily adopted to any observable that 
is a sum of terms with bounded support. To compute 
tr(A|*)(*|) = Ea<btr(AbPab) + Eatr(Aap,) it is Suffi- 
cient to determine the reduced density operators pab and 

Pa- 



For a single WGS (m = l)we obtain pi2 = {Ui (E) 
U2){j:rsM^){t\)iUi(g>U2y with 

N 

rs,t = e-*^ n + e'^=+'^=*"^^-i(^=-*=)i^") (5) 

c=3 

and 7 = J2l,b=l ^ab{SaSb - tah) + Y.l=l{daSa + d'ata). 

This generalizes the formula for WGS without deforma- 
tion obtained in Ref. P|. Eq. demonstrates that for 
any WGS, the reduced density operator of two (and one) 
spins can be calculated with a number of operations that 
is linear in the system size N, as opposed to an exponen- 
tial cost for a general state. 

A straight-forward generalization of Eq. (|SJ) allows one 
to calculate two-qubit reduced density matrices for su- 
perpositions of the form of Eq. Q in time 0{m?N). 
Therefore the expectation value of an observable A of the 
form of Eq. with K terms requires 0{m? KN) steps. 
This implies that even for Hamiltonians where all spins 
interact pairwise (and randomly), i.e. K = N{N — l)/2, 
the expectation value of the energy for our ansatz states 
can be obtained in 0{m?N^) steps. For short-range in- 
teraction Hamiltonians, this reduces to 0{m?N'^). The 
total number of parameters (and memory cost) scales as 
0{N'^+mN), which can be further reduced by employing 
symmetries. 

Regarding (iii) and (iv), one can easily adopt a WGS 
to any given geometry by a proper choice/restriction of 
the adjacency matrix F. A state corresponding to a linear 
cluster state ^3\, for instance, will have only Fa.a+i 0, 
while Ta^a+i would correspond to longer-ranged cor- 
relations. Different choices of (pab lead to very different 
(entanglement) properties: For (pab = |xa — Xf,|~'', where 
Xa denotes the spatial coordinates of spin a, one obtains 
diverging correlation length for two-point correlations, 
while block-wise entanglement can either be bounded or 
grow unboundedly, depending on the choice of /? Q . Sim- 
ilarly, more complicated geometries on lattices in higher 
spatial dimensions can be chosen. 

Variational method. Any state of the form Eq. 
with TO = poly(A^) permits the efficient calculation of 
expectation values of any two-body Hamiltonian H. A 
good approximation to the ground state is then obtained 
by numerical optimization of the parameters character- 
izing the state. Starting from random parameters, one 
descends to the nearest minimum using a general local 
mmimizer (we used L-BFGS [13). Another approach 
that we found to work well is to keep all parameters fixed 
except for cither those corresponding to (i) one local uni- 
tary Ua, (ii) one phase gate Uab{fab) or (iii) the defor- 
mation vector da for one site a. In each case, the energy 
as a function of this subset of parameters turns out to be 
a quotient of quadratic forms, which can be optimized 
using the generalized-eigenvalue (Rayleigh) method. A 
similar result holds for the superposition coefficients aj. 
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One then optimizes with respect to these subsets of pa- 
rameters in turns until convergence is achieved. If one in- 
creases TO stepwise, one — somewhat surprisingly — does 
not get stuck in local minima. 

A significant reduction of the number of parameters 
and the computational costs may be achieved by exploit- 
ing symmetries, or by adapting F to reflect the geometri- 
cal situation. For instance, for systems with short range 
interactions and finite correlation length, one might re- 
strict the range of the weighted graph, i.e. Fab = if 
|xa — Xb| > tq. This reduces the number of parameters 
describing the WGS from OiN"^) to 0{N). For transla- 
tionally invariant Hamiltonians, a better scheme is to let 
Tab depend only on |xa — x^j. This reduces the number 
of parameters to 0{N) as well, and it seems to hardly 
affect the accuracy of the ground state approximation. 
Hence, it allows one to reach high numbers of spins N 
and thus to study also 2D and 3D systems of significant 
size. Trading accuracy for high speed one may even use 
a fully translation-invariant ansatz, where also Da and 
Ua are constant and independent of a. In the latter case, 
for Hamiltonians with only nearest-neighbor interactions, 
the expectation value of the energy can be obtained by 
calculating only a single reduced density operator, and 
the computational cost to treat 2D [and 3D] systems of 
size N = L'^ [N = L^] turns out to be of 0(L) rather 
than 0{N). 

Demonstration. The Ising model. Our method allows 
us to determine, with only moderate computational cost, 
an upper bound on the ground state energy of a strongly 
interacting system of arbitrary geometry. Together with 
the Anderson lower bound, one can hence obtain a quite 
narrow interval for the ground state energy and observe 
qualitative features of the ground state [llj . To illustrate 
our method, we have applied it to the Ising model in ID, 
2D and 3D with periodic boundary conditions, described 
by the Hamiltonian 



H 



(a.b) 



(6) 



where (a, b) denotes nearest neighbors. For a spin chain 
with = 20, and a 2D lattice of size 4 x 4 we compared 
our numerical ground state approximation with exact re- 
sults (Fig.nji). We have also performed calculations for 
larger 2D systems up to 14 x 14. We note that the accu- 
racy can be further improved by increasing to (see Fig. 
^p) . In fact our numerical results suggest an exponential 
improvement with to. We have also tested the fully trans- 
lation invariant ansatz with distance dependent phases, 
constant da and alternating Ua for ID, 2D and 3D sys- 
tems of size A^ = 30, AT = 900 and A^ ^ 27000 respec- 
tively (see Fig. There, for lack of a refernce value 
for the exact ground state, we compare with the An- 
derson bound obtained by calculating the exact ground 
state energy Ea for system size A' = 15,3^,2'^ respec- 
tively. In the 2D and especially the 3D case it is not 



expected that the Anderson bound is particularly tight 
and may lead to a significantly underestimation of the 
precisions achieved by our approach. The states approx- 
imated with this simple ansatz also show qualitatively 
essential features of the exact ground state. As an exam- 
ple, the maximal two-point correlation function QJJja^^ 
(where the two point correlation functions are defined as 

Q'a^p = i^a^'^p'^) - is plottcd agaiust the 

magnetic field B in Fig.[5l3. Strong indication for the oc- 
currence of a phase transition can be observed: the cor- 
relations significantly increase around B « 1.1, 3.12, 5.22 
in ID, 2D, 3D respectively. This is in good agreement 
with estimates employing sophisticated power series ex- 
pansions for the infinite systems or Fade approximants 
based on large scale numerical simulations, which expect 
the critical points at B = 1,3.04,5.14 We also 

remark that the approximated states show a scaling of 
block-wise entanglement proportional to the surface of 
the block, i.e. S'at^ « /?_bL'^"°~^, where (3b is some con- 
stant depending on magnetic field B, Na = L^™ and 
dim is the spatial dimension. We can estimate Pb and 
find that it significantly increases near the critical point. 
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FIG. 1: (Color online.) (a) Relative deviation from exact 
ground state energy for Ising chain with A = 20 (blue) and 
4x4 2D lattice (green) with periodic boundary conditions 
as function of magnetic field B (calculated using BFGS min- 
imization with symmetrized phases, m < 6). (b) ID Ising 
chain with A = 20. Improvement of relative deviation from 
ground state energy as function of number of superposed 
states m for various field values B (calculated using Rayleigh 
minimization without symmetrized phases). 



Generalizations: Our approach can be adapted di- 
rectly to spin-^ systems using the representation Eq. ||2Jl. 
There the sum over binary vectors s with Si = 0, 1 has 
to be changed to n-ary vectors s with = 0, 1, n — 1 
and the corresponding matrices/vectors F, d, J7 have to 
be modified accordingly. However, the limit n — > cx) to 
infinite dimensional systems is both problematic and im- 
practical, as the computational effort increases with n. 
For continuous variable systems we thus choose a closely 
related but slightly different approach. 

The description of field theories on lattices generally 
leads to infinite-dimensional subsystems such as har- 
monic oscillators. A Klein-Gordon field on a lattice for 
example possesses a Hamiltonian quadratic in position 
and momentum operators X and P whose ground state 
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FIG. 2: (Color online.) Ising model in ID (blue) with iV = 30, 
2D (green) with = 30 x 30 = 900 and 3D (red) with 
= 30 X 30 X 30 = 27000 spins arranged as chain, square, and 
cubic lattice, respectively, for fully symmetric ansatz states 
with ipab ~ .f{\xa — a;;, I), da — 1 as function of magnetic field 
B/dim, where dim is dimension of lattice, (a) Relative devia- 
tion of ground state energy {Emf — E)/Emf per bond from 
to mean field approximation Emf (solid), and of Anderson 
bound {Emf — Ea)/Emf (dashed). Translational invariance 
is reduced by using Ui 7^ U2 (alternating), (b) maximal two- 



point correlation 



for nearest neighbors. 



is Gaussian This suggests that techniques from the 
theory of Gaussian state entanglement (see |0| for more 
details) provide the most natural setting for these prob- 
lems. To this end, consider N harmonic oscillators and 
the vector, R= (i?i, i?2Ar)^ = {Xi, Pi, Xn , PnV ■ 
The canonical commutation relations then take the form 
[Rj,Rk] — ifTjk with the symplectic matrix a. All 
information contained in a quantum state p can then 
be expressed equivalently in terms of the characteris- 
tic function XpU) = tr[pW^(^)] where $, G R^^ and 
^(^) = exp(z^-^(TR). Then, expectation values of poly- 
nomials of X and P can be obtained as derivatives of 
X- For Gaussian states, i.e. states whose characteris- 
tic function is a Gaussian Xp(^) = Xp(0)e~*^ €^ 
where here 7 is a 2N x 2Af-matrix and D G is a 

vector, these expectation values can be expressed effi- 
ciently as polynomials in 7 and D. On the level of wave 
functions a pure Gaussian state is given by \F, G; a) = 
ix (F-(G)x+a x|^\ .^^j^gj-g p and G are real 



symmetric matrices, a is a vector, C is the normalization 
and 



F 



-GF-^G -GF-'^ 



D 



GF-ia 



(7) 



Now, we may consider coherent superpositions ~ 
^^ai\Gi, Fi]a.i) to obtain refined approximations of a 
ground state. These do not possess a Gaussian character- 
istic function but a lengthy yet straightforward compu- 
tation reveals that the corresponding characteristic func- 
tion X|V')(V'|(^) ^ smn. of Gaussian functions with com- 
plex weights. Then it is immediately evident that in this 
description we retain the ability of efficient evaluation of 
all expectation values of polynomials in X and P. This 
allows one to establish an efficient algorithm for the ap- 
proximation of ground state properties of lattice Hamil- 



tonians that are polynomial in X and P. 

Summary and Outlook: We have introduced a new 
variational method based on deformed weighted graph 
states to determine approximations to ground states of 
strongly interacting spin systems. The possibility to com- 
pute expectation values of local observables efficiently, to- 
gether with entanglement features similar to those found 
in critical systems, make these states promising candi- 
dates to approximate essential features of ground states 
for systems with short range interactions in arbitrary ge- 
ometries and spatial dimensions. One can also generalize 
this approach to describe the dynamics of such systems, 
systems with long range interactions, disordered systems, 
dissipative systems, systems at finite temperature and 
with infinite dimensional constituents. In fact, general- 
izations of our method that deal with these issues are 
possible and will be reported elsewhere. 
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